{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# *WIP:* Using the `vista` Package\n",
    "\n",
    "This notebook is a work in progress to demo how PVGeo can be used with `vista` for creating integrated visualizations directly in a Python environment. At this time, the 3D rendering is perfromed in a separate window and we have yet to embed the VTK rendering windowinf into some sort of Jupyter widget.\n",
    "\n",
    "Maybe someone reading this knows how to embed the rendering window into a widget?!?? If so, join us on [slack](http://slack.pvgeo.org) and let's collaborate!\n",
    "\n",
    "**DISCLAIMER:** This currently only works on the latest version of `vista` and `PVGeo`\n",
    "\n",
    "```sh\n",
    "pip install -U vista\n",
    "```\n",
    "\n",
    "**Note:** *If not on a webserver*, you can have interactive 3D renderings by specifying `notebook=False` in the plotting routine. This will open a seprate window with the rendering.\n",
    "\n",
    "## Overview\n",
    "\n",
    "The goal set forth for this notebook is to use a comination of Python packages to create an integrated visualization of some data and models for a specific project. These packages and their tasks are:\n",
    "\n",
    "- `discretize` for some file IO,\n",
    "- `SimPEG` inversion results for an inverted model\n",
    "- `PVGeo` for its post processing filters and data integration algorithms\n",
    "- `vista` to create the 3D renderings of the whole data scene\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import vista\n",
    "import PVGeo\n",
    "import numpy as np\n",
    "import discretize\n",
    "\n",
    "print('NumPy Version: %s' % np.__version__)\n",
    "print('PVGeo Version: %s' % PVGeo.__version__)\n",
    "print('vista Version: %s' % vista.__version__)\n",
    "print('discretize Version: %s' % discretize.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import some specific algorithms from PVGeo that we'd like to use\n",
    "from PVGeo.grids import ExtractTopography\n",
    "from PVGeo.ubc import TopoReader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This sets the plotting theme of `vista`\n",
    "vista.set_plot_theme('document')\n",
    "vista.rcParams['cmap'] = 'bwr_r'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the Data\n",
    "\n",
    "Here we load in some data we'd like to process: the [**Laguna del Maule Bouguer Gravity**](http://docs.simpeg.xyz/content/examples/04-grav/plot_laguna_del_maule_inversion.html#sphx-glr-content-examples-04-grav-plot-laguna-del-maule-inversion-py) example from the SimPEG docs.\n",
    "\n",
    "This data scene was produced from the [Laguna del Maule Bouguer Gravity](http://docs.simpeg.xyz/content/examples/04-grav/plot_laguna_del_maule_inversion.html#sphx-glr-content-examples-04-grav-plot-laguna-del-maule-inversion-py) example provided by [Craig Miller](https://github.com/craigmillernz) (Maule volcanic field, Chile. Refer to Miller et al 2016 EPSL for full details.)\n",
    "\n",
    "> Miller, C. A., Williams-Jones, G., Fournier, D., & Witter, J. (2017). 3D gravity inversion and thermodynamic modelling reveal properties of shallow silicic magma reservoir beneath Laguna del Maule, Chile. Earth and Planetary Science Letters, 459, 14–27. https://doi.org/10.1016/j.epsl.2016.11.007"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the a TensorMesh and some already processed model data \n",
    "mesh = discretize.TensorMesh.readUBC('craig_chile.msh', directory='data/Craig-Chile')\n",
    "models = {'lpout': mesh.readModelUBC(fileName='Lpout.mod', directory='data/Craig-Chile')}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load topography data using PVGeo\n",
    "topo = TopoReader().Apply('data/Craig-Chile/LdM_topo.topo')\n",
    "# Note that PVGeo will return a vista wrapped data object if vista is available!\n",
    "topo.set_active_scalar('Elevation')\n",
    "topo.plot(cmap='gist_earth', clim=[1.7e+03, 3.104e+03])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build a Pipeline \n",
    "\n",
    "Here we build up a pipeline that will transform and integrate the data. This pipeline uses several algorithms taht pass their output onto the next algorithm and so forth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This adds an active cell array to describe how the topo surface splits the dataset\n",
    "extractor = ExtractTopography(offset=-150, tolerance=10, op='underneath')\n",
    "extracted = extractor.Apply(mesh.toVTK(models=models), topo)\n",
    "extracted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `ExtractTopography` algorithm in the above cell adds a new array called `'Extracted'` that has 0s and 1s for how the topography surface splits the dataset. We can now use this array to extract the region of the model that is beneath the topography (the subsurface region of the model).\n",
    "\n",
    "To do this, we'll use a simple `PercentThreshold` filter from PVGeo. The `PercentThreshold` algorithm defaults to threshold at 50% so it will remove all the cells with 0s and preserve the region of the model where the cells are active (the 1s). There are many ways you could threshold this model to extract the subsurface but this algorithm provides a clean way to do it with default parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# threshold out the topography by a simple percent threshold in PVGeo\n",
    "subsurface = extracted.threshold_percent(scalars='Extracted')\n",
    "subsurface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the output of the above cell now shows a `vtkUnstructuredGrid`. This dataset contains only the cells that were extracted by the `ExtractTopography` algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the scalar range to plot everything the same way\n",
    "rng = subsurface.get_data_range('lpout')\n",
    "# this range helps us create consistent and meaningful color legends"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the subsurface region of the model domain, we can use that dataset for further volume extraction and slicing. The next cell uses a filter directly from `vista` that will threshold a specific value range of the dataset so that we will have a volumetric body to inspect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# threshold using vista on a specific value\n",
    "low = subsurface.threshold_percent(35, invert=True)\n",
    "low.plot(rng=rng, window_size=[1024//2, 768//2])\n",
    "\n",
    "high = subsurface.threshold_percent(85)\n",
    "high.plot(rng=rng, window_size=[1024//2, 768//2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Or better yet, use an interactive tool\n",
    "#vista.Threshold(subsurface)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the `ManySlicesAlongAxis` filter from PVGeo to create `n` slices of the subsurface model dataset along a given axis. In our case, we'll slice along the Y-axis (index 1) and create 5 slices.\n",
    "\n",
    "Note that we can interactively plot any of these datasets by ensuring the `notebook` argument is set to `False` which will open a seperate rendering window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Slice the subsurface model\n",
    "slices = subsurface.slice_along_axis(5)\n",
    "slices.plot(rng=rng, notebook=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot in an Interactive 3D Window\n",
    "\n",
    "Okay, so now we have all these datasets and filter products of those datasets... what if we need to spatially relat all of those datasets and interactively view them? For example, lets display the thresholded volume with the slices of the data set all underneath the topography surface:\n",
    "\n",
    "\n",
    "First, we create the plotter/rendering environment:\n",
    "\n",
    "```py\n",
    "p = vista.Plotter(notebook=False)\n",
    "```\n",
    "\n",
    "Next we add all the datasets we'd like to display by repeatedly calling `p.add_mesh` with the data object and any plotting arguments such as colormap or the scalar array you'd like to color it by:\n",
    "\n",
    "```py\n",
    "cmap = 'bwr_r'\n",
    "p.add_mesh(low, scalars='lpout', rng=rng, showedges=False)\n",
    "p.add_mesh(high, scalars='lpout', rng=rng, showedges=False)\n",
    "p.add_mesh(slices, scalars='lpout', rng=rng, showedges=False, opacity=0.85)\n",
    "p.add_mesh(topo, opacity=0.5, psize=1.0, colormap='gist_earth', rng=[1.7e+03, 3.104e+03])\n",
    "```\n",
    "\n",
    "Then we add any extra plotting labels such as the axis orientation and spatial labels/bounds:\n",
    "\n",
    "```py\n",
    "p.add_bounds_axes(extracted, color='k', fontsize=30)\n",
    "p.add_axes()\n",
    "```\n",
    "\n",
    "Lastly, we call `p.show` to display the rendering windoe much like how you might use `plt.show` in `matplotlib`.\n",
    "\n",
    "```py\n",
    "p.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot them all in one rendering environment!\n",
    "p = vista.Plotter()\n",
    "\n",
    "#p.add_mesh(extracted, syle='wireframe', opacity=0.25, scalars='lpout', clim=rng)\n",
    "p.add_mesh(low, scalars='lpout', clim=rng)\n",
    "p.add_mesh(high, scalars='lpout', clim=rng)\n",
    "p.add_mesh(slices, scalars='lpout', clim=rng, opacity=0.85)\n",
    "p.add_mesh(topo, opacity=0.5, point_size=1.0, cmap='gist_earth', clim=[1.7e+03, 3.104e+03])\n",
    "\n",
    "p.show_bounds(color='k', font_size=30, location='outer')\n",
    "p.add_axes()\n",
    "\n",
    "p.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:dev]",
   "language": "python",
   "name": "dev"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
